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The mutual neutralization of and H~ at low collision energies is studied by means of a 
molecular close-coupling approach. All degrees of freedom are treated at the full quantum level also 
taking into account the identity of the nuclei. The relevant ^E^ and electronic states as well as 
the associated non-adiabatic radial couplings are calculated for internuclear distances between 0.5 
and 50 a^. Following a transformation into a strictly diabatic basis, these quantities enter into a 
set of coupled equations for the motion of the nuclei. Numerical solution of these equations allows 
the cross sections for neutralization into the H(l) + H(n), n = 1,2,3 final states to be calculated. 
In the present paper, results are reported for the collision energy region 0.001 to 100 eV, with 
special emphasis on the important energy region below 10 eV. The low temperature rate coefficient 
is obtained from a parametrization of the calculated cross section and is estimated to be valid over 
the range 10 to 10,000 K. 

PACS numbers: 34.70.+O, 31.50.Df, 31.50.Gh 



I. INTRODUCTION 

The mutual neutralization of and H~ is the pro- 
totype reaction for electron transfer between two oppo- 
sitely charged ions. With only two protons and two elec- 
trons involved, it is the simplest ion-ion reaction possi- 
ble. Here, calculations should be able to provide accurate 
results which could serve as benchmarks in the develop- 
ment of new theoretical methods and models. Further- 
more, the reactants have well defined quantum states and 
are relatively easy to prepare in the laboratory. A de- 
tailed comparison between experiment and theory should 
thus be possible and should provide a natural starting 
point to understanding more complicated ion-ion colli- 
sion phenomena. However, the simplicity of this reaction 
is deceptive, a fact which has stimulated a large num- 
ber of studies in the past. Nevertheless, surprisingly few 
of these have investigated the neutralization cross sec- 
tion for collision energies below a few eV, which plays 
an important role in determining the H~ abundance in 
low temperature ionized environments. A compelling ex- 
ample from astrophysics is its influence on the gas phase 
formation of H2, in particular in the chemistry of the 
primordial gas, where one of the main reaction paths in- 
volves H~ as an intermediate [J, H, 

Bates and Lewis ^4^ were the first to study the mutual 
neutralization of H^ and H^ , which is written schemat- 
ically as 



H^ 



H" 



H(l) +H(n) 
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where n is the main quantum number of the excited hy- 
drogen atom. Their calculations were based on the semi- 
classical Landau-Zener model @| in which the reaction 
is thought to proceed along a system of horizontal co- 
valent potential energy curves crossed by an attractive 
ion-pair potential. This pioneering study was later fol- 
lowed by others, among them by Olson et al. [3] and more 
recently by Eerden et al. [§] . Fussen and Kubach have 
gone beyond the Landau-Zener concept and investigated 
the reaction at low collision energies using a one electron 
close-coupling model. 

The only measurement of the neutralization cross- 
section for collision energies below 3 eV is the merged 
beam experiment of Moseley et al. [l^]- Their results 
seem to indicate a cross section that is approximately 
a factor of three higher than obtained in most of the 
cited theoretical studies. However, the theoretical results 
are mutually consistent and, for higher energies, also in 
more or less good agreement with the measurements re- 
ported by Szucs et al. [lH and Peart and Hayton p^ . 
This intriguing situation has led to the suspicion that 
the low en ergy cross section of Moseley et al. is overesti- 
mated 

In this paper we report results of a fully quantum me- 
chanical ah initio study of the reaction at hand. We are 
mainly concerned with the collision energy region 0.001 
to 10 eV, but the calculations have been extended up to 
100 eV in order to facilitate comparison with the cited 
experiments. The present study is based on a molecu- 
lar close-coupling expansion of the total wave function. 
Accurate electronic structure methods have been used to 
calculate adiabatic states and non-adiabatic couplings of 
both gerade and ungerade symmetry over a wide range 
of internuclear distances. Subsequently, a transforma- 
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tion into a strictly diabatic representation has been per- 
formed and the diabatic quantities entered into a set of 
coupled equations for the motion of the nuclei. By nu- 
merically solving these equations, we have been able to 
calculate the neutralization cross sections for scattering 
within both the gerade and ungerade inversion symme- 
tries. Finally, the identity of the nuclei has been ac- 
counted for by properly combining the gerade and unger- 
ade cross sections. 

The paper is organized as follows. Section HIl describes 
the theoretical framework underlying the present study. 
Section Hill describes the electronic structure calculations 
and discusses the resulting potential energy curves and 
coupling matrix elements. Section IIVI gives a brief ac- 
count of the numerical treatment of the scattering prob- 
lem. SectionlVlpresents and discusses the calculated neu- 
tralization cross sections and the reaction rate coefficient. 
Section IVTl concludes the paper and summarizes the main 
results. Unless stated otherwise, vector and matrix quan- 
tities are written in boldface and atomic units are used. 



II. THEORY 



where E is the total energy in the center of mass system. 
We expand ipiR, r) in terms of electronic and nuclear 
states, 0f(R, r) and xf(R-), according to 



V'(R,r) = 



KR,r)xr(R) = 0"x'' 



(4) 



where (R, r) are normalized solutions to the electronic 
Schrodinger equation 



i/°V-(R,r)=e^(i?)0^(R,r) 



(5) 



The electronic energies ef{R) defend parametrically on R 
and obey the non-crossing rule . The particular choice 
of basis functions defined by Eq. (O is called the adia- 
batic representation and is denoted by the superscript a. 
Inserting the expansion (HI into Eq. ([3]) and using the or- 
thonormality of the electronic states yields a set of cou- 
pled differential equations for the motion of the nuclei: 



2n 



R 



.pa . y. 



R 



2/z 



(6) 



The molecular expansion: adiabatic and 
diabatic formulations 



Consider a diatomic system composed of two nuclei 
and N electrons. The masses and charges of the nuclei 
are denoted by Mi and Zi, respectively. Let R = (i?, 9, ip) 
be the internuclear separation vector and r.^ a vector 
pointing from the center of mass of the nuclei to the 
z-th electron. All vectors are expressed with respect to 
a space fixed axis system. For convenience, the notation 
r = (ri,r2, ■ • ■ ,J^n) will occasionally be used. In these 
coordinates, the center of mass motion of the total sys- 
tem can be separated out, and the Hamiltonian takes the 
form UM 



with 



2n ^ 



(1) 



where /i is the reduced mass of the nuclei and is 
the electronic Hamiltonian. Neglecting terms of the form 
[2 (Ml + M2)]"^Vri • Vr^., it follows that 



2 ' R 



N 

E 

i>j=i 



r, - r,- 



N 

E 



Z2 



h^R + r, 



(2) 



The total wave function iJ^CR, r) satisfies the time- 
independent Schrodinger equation 



H'iJj{R,r) = E'4j(R,r) 



(3) 



Ft^.(R) = ((/)nVR|0|) 



Vt,(R)^{cb^\H^'\cl>^) 



(7) 



(8) 



(9) 



By definition, the matrix is diagonal and contains 
the adiabatic potential energy curves ef(i?). The matrix 
vector and the matrix both have off-diagonal el- 
ements, which are the non-adiabatic couplings between 
different electronic states. In addition, C also has diag- 
onal elements commonly referred to as adiabatic energy 
corrections. If the electronic wave functions are taken to 
be real, the chain rule for derivatives can be used to show 
that is antisymmetric, i.e.. 



F-(R) = -F^,(R) . 



(10) 



The electronic states approximately correlate with 
the separated atomic limits and therefore each term in 
Eq. (HI) defines a scattering channel. To determine the 
outcome of any collision process, we need to solve the 
coupled equations out to the asymptotic region where 
these limits are reached. However, it has been known for 
a long time that a product of the form (R)(/)f (R, r) 
does not properly describe the translational motion of 
the electrons along with the nuclei [13, HI, HI]. This 
can lead to difficulties of both formal and practical na- 
ture. One of the most prominent is the appearance of 
non- vanishing cou pling s at large or even infinite internu- 
clear distances [19, Il9,^0l- However, it is demonstrated 
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in Sec. IV Bl that the influence of these coupUngs on the 
present results is expected to be small. 

The non-adiabatic coupling terms in Eq. ([6]) often vary 
rapidly with R and may cause problems in the numeri- 
cal solution procedure. Furthermore, the presence of F'' 
leads to differential equations which contain both first 
and second order derivatives. It is therefore desirable to 
transform Eq. ^ into an electronic basis where some or 
all of the non-adiabatic couplings disappear. Following 
Mead and Truhlar [2l| , we shall refer to a representation 
in which all three components of F'' vanish as strictly 
diabatic. From this point on, it is assumed that the to- 
tal wave function can be represented in a finite set of M 
adiabatic basis functions. The transformation may then 
be written in terms of an orthogonal matrix Ty(R) 
as 



M 



(R,r)7;^| 



R) 



(11) 



In order for the total wave function to be preserved, the 
nuclear states must transform according to 



M 

xf(R)=5]x|(R)?;T(R) • 



(12) 



Using the relation (I12p to replace the adiabatic states in 
Eq. ^ yields after some manipulation 



1 

2^" 



1 



1 



IV^ - -F'^ • Vr - -^G^ + U'' 
H 2fi J 

- Elx'' , (13) 
F'^ = T"^ (Vr + F^) T , (14) 



with 



2F'' • V 



R 



G^IT 



(15) 



(16) 



The first derivative coupling F'^ expressed in the new 
basis is seen to vanish provided that T is a solution to 
the equation 



(Vr + F^)T = 



(17) 



In the special case where all electronic states have the 
same angular momenta projection onto the molecular 
axis, only the radial component of F^ can be non- 
zero [21]. Equation p7p then reduces to 



where 



d 

dR 



T(i?) 



T(i?) = 



dR 



(18) 



(19) 



Assuming that the first derivative couplings between the 
adiabatic states with i<M and those with i > M are 
zero, also G'^ can be shown to vanish [22] • The coupled 
equations then take the form 



1 9 



2fi 



(20) 



In the strictly diabatic representation, the non- 
adiabatic couplings transform into off-diagonal elements 
of the potential matrix (electronic couplings) and as 
a consequence the diagonal elements of this matrix (di- 
abatic potentials) no longer obey the non-crossing rule. 
As will be illustrated in Sec. IIII CI the strictly diabatic 
representation may show little resemblance to the simple 
and intuitive curve crossing picture often referred to as 
diabatic. 



B. Asymptotic boundary conditions and the 
integral cross section 

The interpretation of xf(R) is usually simplified by 
requiring that 



hm T(i?) 

it— >00 



(21) 



which makes the adiabatic and the diabatic representa- 
tions coincide in the region where the scattering wave 
function is evaluated. Thus, in any of the representa- 
tions, the physical situation is that of an incoming plane 
wave in some entrance channel j, and outgoing spheri- 
cal waves in all channels which are energetically allowed. 
The direction of the incoming plane wave can be chosen 
along the space fixed z-axis. If the potential falls off rea- 
sonably fast (23j . the appropriate boundary conditions 
are given by 



xf(R)^^ S.,,e 



(22) 



where {E, 9) is the scattering amplitude for entering 
the collision in channel j and leaving it in channel i. Since 
the potential matrix is spherically symmetric, the scat- 
tering amplitude is independent of the azimuthal angle ip. 
The asymptotic wave number ki is defined as 



,1/2 



(23) 



where Ef" = limi^.^oo VuiR) is the corresponding thresh- 
old energy. In the case of a Coulomb potential, the ex- 
ponents in Eq. (|22p have to be modified with logarithmic 
phase factors which arise due to the long range nature of 
this interaction [23]. 

Evaluating the probability fluxes associated with the 
two terms in Eq. leads to the familiar multichannel 
expression for the integral cross section: 



a,,{E) 



27rfc, 



r\ME, 

Jo 



■sinede 



(24) 
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If the adiabatic states are connected by non- vanishing 
couphngs at infinity, neither the boundary condition (PT|) 
nor the scattering formahsm developed below will be ap- 
propriate. In this case, the couplings can be manually 
cut at some large but finite internuclear distance. This 
procedure is usually justified if the collision energy is low 
enough, but its validity should be considered from case 
to case. 



for the scattering amplitude in terms of the open-open 
or physical partition of S; . The integral cross section for 
scattering from channel j to channel i is readily evaluated 
to be 



where 



(31) 



C. The scattering formalism 



(32) 



The conventional approach to solving Eq. ([20|l is to 
express x'^ in terms of a partial wave expansion. 



oo 

xf(R)--^A;K,/i?)PKc 



(25) 



where Pi{cos9) are the well known Legendre polynomi- 
als [13] and Ai are constants to be chosen so that the 
boundary conditions l|22p are fulfilled. The radial wave 
functions Ui^R) can be shown to satisfy the equations 



1 d" 



2^di?2 2^i?2 



M 



«.,,(i?) + 5]l/,^U,,(^) 



Eu^^R) , (26) 



subject to the physical boundary conditions Ui.;(0) = 0. 
If the corresponding derivatives u[ ;(0) are chosen appro- 
priately, these equations have M linearly independent 
solutions which can be combined into a square matrix u. 
It follows that 



di?2- 



1 + Qz = , 



where 



Qi = 2^l {El - V^) - 



^(^ + 1) 

i?2 



(27) 



(28) 



It is now assumed that the potential matrix is diagonal 
beyond some radius p and here has elements correspond- 
ing to either pure Coulomb or pure short-range interac- 
tions. In this region, Eq. (|27p is analytically solvable and 
the radial wave functions may be expressed in terms of 
the incoming and outgoing wave solutions, aij^i{R) and 
(3ij,i{R), as defined in the Appendix. These are essen- 
tially the same as adopted by Johnson [2^ but have been 
extended to cover the Coulomb case. Thus, 



ui{R) = aiiR) - f3iiR)Si 



(29) 



which defines the scattering matrix S/ . Inserting this ex- 
pression into the partial wave expansion (|25p and com- 
paring with the boundary conditions (p2)) leads to the 
well known formula 



1 



2i(fc,fcj)^/2 



1=0 



(30) 



Instead of working with S; directly, it is convenient to 
introduce the reactance matrix K; [26j which is defined 
by 



uiiR) - a,(i?)-bj(i?)K; 



R>p 



(33) 



where aijj{R) and bij^R) are regular and irregular solu- 
tions to Eq. ((27)) explicitly given in the Appendix. Com- 
bining the definitions above together with those from the 
Appendix leads to the relationship [25] 



sr = (i+iKn 



zKr) 



(34) 



between the open-open partitions of S; and K;. This ex- 
pression resembles the Cayley transform and, since S?° is 
both symmetric and unitary, implies that K?° is real [2 



Finally, we introduce the logarithmic derivative 



(35) 



which transforms the radial Schrodinger equation into 
the numerically more stable matrix Riccati equation (27l . 



yj + Q/ + y? = , 



(36) 



while the boundary condition at the origin becomes that 
of a diagonal matrix with infinite elements. This is the 
set of coupled equations that we ultimately have to solve. 
The relationship between y; and K; is easily seen to be 



K, - [yi{p)hi{p)~h[{p)] 
X [yi(p)a/(p) -a;(p)] 



(37) 



D. Modifications in the case of identical nuclei 

Up to this point in the discussion complications that 
might arise if the nuclei are identical have been ignored. 
These may be considered in two steps. Firstly, since the 
charges of the nuclei that enter H'^^ are equal, the elec- 
tronic states will be either symmetric (gerade) or anti- 
symmetric (ungerade) under inversion of r. If the nuclei 
are labeled A and B, and a,/?, 7, ... represent groups of 
electrons in some definite quantum states, each asymp- 
totic channel will be a linear combination of two config- 
urations, say A(a) -I- B(/3) and A(/3) -I- B(a). Such linear 
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combinations can not describe the situation where a nu- 
cleus is known to carry a certain number of electrons in a 
certain quantum state (unless a — (3). Localization of the 
electron cloud can be achieved by forming proper linear 
combinations of the gerade and ungerade solutions. The 
g/u symmetry is explicitly indicated as a superscript and 
the channel index i is allowed to include all other quan- 
tum numbers. The appropriate scattering amplitudes are 
then given by 

ffi{E,e) = ^[ff^{E,9) + f-{E,e)] (38) 

and 

f!^iE,9)^^[f^^^{E,e)-fr,iE,e)] , (39) 

corresponding to the direct and exchange reactions 

A(a)+B(/3) ^ A(7) + B((5) (40) 

and 

A(a) B(/3) ^ A(^) + B(7) . (41) 

Here, the amplitudes f^j and J^" are given by the usual 
formula ((30)) applied to the gerade and ungerade mani- 
folds separately. 

A second complication arises from the fact that also 
the masses of the nuclei are equal. Consequently, it is 
impossible to distinguish direct scattering in the direc- 
tion 6 from exchange scattering in the direction n — 6. 
This effect can be accounted for by adding the direct and 
exchange amplitudes coherently. The fully symmetrized 
scattering amplitudes are [1^ 

f^{E, 6) = f^{E, 6) ± /-(i?, ^ - 0) , (42) 

where the sign depends on whether the spatial part of 
the wave function should be symmetric (+) or antisym- 
metric (-) under exchange of the nuclei. To illustrate this 
point we consider a hydrogen quasi molecule where the 
two protons are known to be in a singlet spin state. Ac- 
cording to the Pauli principle, the total wave function 
(spin part included) must be antisymmetric under ex- 
change of the nuclei. Since the singlet spin state itself is 
antisymmetric under such an operation, the spatial part 
of the wave function is forced to be symmetric. Thus, f^j 
is the appropriate choice of amplitude. 

The integral cross section is obtained by taking the 
modulus square of Eq. ([i^ and integrating the result 
over the unit sphere: 

4(^)= E <AE)+Y.^UE) (43) 

I even / odd 

and 

I odd I even 



where crf^- ; and cr^_; are given by Eq. (|32|) applied to the 
gerade and ungerade manifolds separately. The subscript 
even (odd) indicates that only even (odd) partial waves 
should be summed over. In the particular case of the hy- 
drogen quasi molecule, the total spin state of the nuclei 
is either a singlet (para hydrogen) or a triplet (ortho hy- 
drogen). Taking into account the degeneracy factors (1 
and 3) and the exchange symmetry of these spin states, 
the symmetrized integral cross section for an ensemble of 
particles with randomly oriented spins is 

-r(^) - jE4.(^) + iE4..(^) 

I even / odd 

+\ E <^(^) + JE<'(^)- 

I even I odd 

(45) 

This can be compared with the expression 

^t'iE) = \ Y.<AE) + \Y.''UE)^ (46) 
/=o /=o 

obtained by treating the nuclei as distinguishable, and 
thus adding the direct and exchange amplitudes incoher- 
ently. 

In deriving Eqs. ([^5]) and (HH]) we have assumed that all 
final states exist in both gerade and ungerade versions. 
This is not always the case. Of particular interest to the 
present study are processes in which the final electronic 
state corresponds to 7 = 5 and exists only in a gerade 
version. In this case, Eqs. (|45)) and (|46)) are still valid 
(setting (Tj" ; = 0) , although the arguments leading to 
these formulas are somewhat modified. 



III. ELECTRONIC STRUCTURE 
CALCULATIONS 

The electronic states relevant to the present study are 
those of H2 . The amount of data published for this sys- 
tem is of course substantial. Of particular relevance here 
are the very accurate electronic structure calculations re- 
ported by Wolniewicz and co-workers [13, 31, 32, 33113^. 
|35| . [36| using explicitly correlated basis functions and the 
calculations by Detmer et al. [13] using a non-spherical 
Gaussian basis set. Reactions like the present one, how- 
ever, require the calculation of not only a large number 
of excited states but also the associated non-adiabatic 
couplings for a wide range of internuclear distances. To 
the best of our knowledge, a complete and accurate set 
of data for these quantities does not exist in the litera- 
ture. We have therefore performed ab initio electronic 
structure calculations of all the relevant potential energy 
curves and non-adiabatic couplings. Where such data is 
available, our results have been cornpared with the bench- 
mark calculations of Refs. [13, [M Isl, HI [13, [H, HI • 
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TABLE I: Specification of the spherical Gaussian basis set 
used in the present study. 



Type 


Exponent 


Type 


Exponent 


s" 


- 


P 


0.0275700 


s 


0.7977000 


P 


0.0094220 


s 


0.2581000 


P 


0.0031910 


s 


0.0898900 


P 


0.0010470 


s 


0.0251300 


d 


2.0620000 




n 0078770 


d 


n 6620000 

yj .\J\J ^UKJXJKJ 


s 


0.0029256 


d 


0.1900000 


s 


0.0008852 


d 


0.0635300 


s 


0.0002917 


d 


0.0211100 


p 


2.2920000 


d 


0.0070570 


p 


0.8380000 


d 


0.0023460 


p 


0.2920000 


f 


1.3970000 


p 


0.0848000 


f 


0.3600000 



"Contraction of six primitive s functions with exponents (82.64, 
12.41, 2.824, 0.7977, 0.2581, 0.08989) and coefficients (0.002006, 
0.015343, 0.075579, 0.256875, 0.497368, 0.296133). 



A. Potential energy curves 

We have calculated the adiabatic potential energy 
curves and the associated electronic wave functions corre- 
sponding to the seven lowest states and the six low- 
est states, here denoted as (1-7)^I]+ and (1-6)^S+, 
respectively. All calculations have been performed at the 
full configuration interaction (FCI) level using a modi- 
fied version of the Dalton 2.0 program package The 
molecular orbitals have been obtained in a {lls,8p,7d,2f) 
spherical Gaussian basis set contracted to [9s,8p,7d,2f]. 
The compact basis functions of this set have been taken 
from the augmented correlation-consistent polarized va- 
lence quadruple zeta (aug-cc-pVQZ) basis set of Dunning 
and co-workers [s^ \^ , while the more diffuse ones have 
been approximately optimized to obtain a good repre- 
sentation of the excited states. The basis set is given in 
Table H 

Figure [2 shows the (2-7)^1;+ and (1-6)^^,+ potential 
energy curves in the range 0.5 to 50 ap. As is well known, 
only one bound state of the negative hydrogen ion exists; 
the ^S*^ ground state [4l[. Accordingly, there is one ger- 
ade and one ungerade electronic state which asymptoti- 
cally correlate with the H"'' -I- H~ ion-pair limit. As the 
internuclear distance is decreased from infinity, the ion- 
pair configuration is successively carried on through a se- 
ries of avoided crossings involving states correlated with 
the H(l) -I- H(ri) covalent limits. In the neighborhood of 
these crossings, the first derivative radial couplings (fT51) 
are significant and are able to cause transitions between 
different adiabatic states. 

Beginning with the gerade symmetry, the (1-6)^!]^]" 
manifold consists of states correlated with the n — 1,2,3 
covalent limits. Corresponding to each limit H(l) + H{n) 
there are n of these states in total. The 7^E+ state has 
the ion-pair configuration from 36 ag up to approximately 
280 ao- Here, the adiabatic wave function changes char- 
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1 / (b) 



Internuclear distance (units of ao) 

FIG. 1: Adiabatic potential energy curves for (a) the (2-7)^E+ 
states and (b) the (1-6)^E^ states of H2. The separated 
atomic limits are indicated at the far right of each subfig- 
ure. For a discussion of the particular ion-pair limit, see the 
main text. 



acter into the n — \ covalent configuration while the ion- 
pair configuration is carried on to another higher excited 
state. At such large distances, however, the motion of 
the nuclei is highly diabatic and the change of character 
of the adiabatic wave functions will not induce any sig- 
nificant ionic-covalent transitions 0, H^, 13 ■ With these 
considerations in mind, we will refer to the 7^1]^ state 
as that which is correlated with the ion-pair limit. 

In the ungerade symmetry, the (l-5)^Ej manifold con- 
sists of states correlated with the n = 2, 3 covalent limits 
while the G^Ej state, with the same considerations as 
above, can be said to correlate with the ion-pair limit. 
Because of symmetry restrictions, there can be no sin- 
glet ungerade state dissociating into the n — \ covalent 
limit. Note that while the gerade and ungerade states 
behave very differently at small internuclear distances, 
these differences vanish quickly as the separation of the 
nuclei is increased. In fact, at 10 ao the energy difference 
for all excited states is less than 9 x 10~^ i?h (Hartree) 
and at 30 ao it is less than 5 x 10~^ £"11. 

Figure [T] illustrates the two avoided crossings near 12 



7 



and 36 Rq, which are related to a change from ionic to 
n = 2 and 3 covalent character (or vice versa) of the 
adiabatic wave functions, respectively. In what follows 
we will refer to them simply as the n — 2 and 3 curve 
crossings. Figure [1] also shows a rich pattern of avoided 
crossings taking place at smaller intcrnuclear distances. 
This region, however, is to a large extent shielded by 
the centrifugal term in Eq. (pS)) and the influence of the 
resulting radial couplings upon the total neutralization 
cross section is rather small (see Sec. IV Bp . The higher 
excited ^S+^j states that are not considered in the present 
study are only accessible to the system through avoided 
crossings at small internuclear distances [s^]. The error 
introduced by excluding these states is therefore expected 
to be of less significance. In principle, also rotational 
couplings to states of ^11 symmetry should be considered. 
However, due to the presence of strong radial couplings 
in the vicinity of the n — 2 and 3 curve crossings, and 
due to the low collision energies considered in the present 
study, these are assumed to be negligible d, Q . 

To evaluate the quality of our calculated potential 
energy curves, we have compared them with those re- 
ported by Wolniewicz and co-workers. Their calcula- 
tions cover the 1^5]+ state out to 12 ag [H, [l^l, the (2- 
6)iE+ states out to 20 ao HI and the (1-6) states 
out to 150 ao d^. The potential energy curve for the 
4^E+ state has been recalculated and extended out to 
80 ao nil . With the energy scale used in Fig. [U these 
potentials would entirely overlap with ours. A more de- 
tailed comparison shows that the largest energy differ- 
ence is obtained for the inner part {R < 1.5 ao) of the 
l^E^ ground state potential, which is of little impor- 
tance to the present study anyway. Here, our computed 
values are approximately 5.7 x lO^'' to 1.5 x 10"'^ i^h 
higher in energy compared to those in Ref. [13]. For 
all of the excited states, our calculated energies are ap- 
proximately 1.0 X 10"'' to 5.6 X 10"'' Ell higher than in 
Refs. m, m, H^. It is of particular interest to estimate 
the quality of the (4-7)-'!]+ and (3-6)^S+ states near the 
critical n = 3 curve crossing. Concerning the unger- 
ade states, the covalent parts of our potentials are within 
1.6 X lO^'' Eh and the ionic parts within 3.0 x lO^'' E^ to 
those reported in Ref. [s^ . Due to the near degeneracy of 
the gerade and ungerade states, we expect this accuracy 
to be common to both inversion symmetries. The qual- 
ity of the adiabatic states and in particular their energy 
splitting at the n — 3 curve crossing will be discussed 
further in the next section. 



B. Radial couplings 

We have calculated the first derivative radial cou- 
plings between all of the adiabatic states considered 
in the previous section. The accurate evaluation of these 
quantities is not trivial. The magnitude and shape of the 
radial couplings are usually very sensitive to the quality 
of the electronic wave functions and it is essential that 



not only the individual states are good in a variational 
sense, but also that the relative energies of these states 
are well represented. In regions where two or more adia- 
batic states become nearly degenerate, further complex- 
ity is added by the fact that a minor change in the atomic 
basis can cause the potential energy curves to artificially 
pseudo cross, giving rise to dramatic effects in the radial 
couplings. These effects do not change the outcome of 
the dynamics but can obscure both interpretation and 
comparison with other results. 

To evaluate the radial couplings we have implemented 
a three point version of the finite difference method de- 
scribed in Ref. [i^. This method provides a systematic 
way to converge the calculations toward the exact result 
in the particular type of basis and wave function being 
considered. We have examined carefully how the cal- 
culated radial couplings depend on the derivative step 
length A_R, as well as the numerical accuracies, Smo f^nd 
6ci, with which the molecular orbitals and the CI wave 
functions are obtained. Stable and converged results were 
observed when AR = 5 x 10~^ ao, Smo — 10"^'' E^ and 
Sci — 10^^^ Eh- The antisymmetry condition pU|) can 
be used as a simple consistency check of the results. In 
the present case, most of the calculated coupling elements 
satisfied this condition to at least four significant digits. 
For practical reasons, it is still desirable to have a cou- 
pling matrix that is fully antisymmetric. This has been 
accomplished by taking each element to be the geomet- 
rical mean of Tij and —Tji. 

The radial couplings are too many to discuss in full 
detail and so we will focus on those relevant for the im- 
portant n = 2 and 3 curve crossings. The notation {i,j)g 
is used to label the radial coupling between the i-th and 
j-th state (and similarly for the ungerade states). 
Due to symmetry, gerade and ungerade states are not 
coupled to each other. 

In Fig. [2] we show the radial couplings among the (2- 
4)^S]+ and (1-3)^S+ states in the neighborhood of the 
n ~ 2 curve crossing. As can be expected from the 
appearance of the potential energy curves, each of the 
(2,4)g and (1,3)„ couplings show a rather broad peak 
near the crossing point at 12 ao resulting from the ex- 
change of ionic and n — 2 covalent character between 
the adiabatic wave functions. Similar shapes and mag- 
nitudes, although displaced to larger internuclear dis- 
tances, are also exhibited by the (2,3)g and (1,2)„ cou- 
plings, which connect the adiabatic states tending to the 
n = 2 covalent limit. The (3, 4)g and (2, 3)„ couplings 
are, on the other hand, much less pronounced in the 
crossing region. For comparison, also plotted in Fig. are 
the results obtained by Wolniewicz and Dressier [soTlssj , 
which include the radial couplings among all of the (2- 
4)^I]+ states and between the 1 and 2^1]+ states. The 
overall agreement with the present results is very good, 
in particular in the region up to approximately 16 ao. 
For larger internuclear distances some deviations can be 
observed. These are most likely due to a slight unbal- 
ance in our description of the nearly degenerate n = 2 
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10 12 14 16 
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FIG. 2; Radial couplings among (a) the (2-4) states and 
(b) the (1-3)^ states in the neighborhood of the n = 2 
curve crossing. The notation {i,j)g label the radial coupling 
between the i-th and j-th state (and similarly for the 
ungerade states). The present results (solid lines) are com- 
pared with the results obtained by Wolniewicz and Dressier 
(crosses) [30l. |33|]. 



adiabatic states. However, as mentioned at the very be- 
ginning of this section, the physical relevance of these 
types of effects is expected to be small. 

Figure[3]illustrates the radial couplings among the adi- 
abatic states associated with the n = 3 curve crossing. 
Here, the gerade and ungerade states are almost com- 
pletely degenerate and so only the results for the for- 
mer of these, i.e., the (4-7) manifold of states, are 
shown. All the couplings appear as peak like structures 
with their maxima located in the region 35.5 to 37.0 aq. 
The most prominent of these is the (4, 7)g coupling which 
reaches its maximum value at 35.9 ap. This coupling 
arises from the exchange of ionic and n = 3 covalent 
character between the 4 and 7^S+ adiabatic states. We 
note that the location of its maximum agrees well with 
the crossing point of the pure ionic and covalent ener gies 
at 36.0 ao obtained from the ground state energy [46| 
and the polarizability [l^l of the H~ ion. The calcula- 
tions of Wolniewicz and Dressier do not extend out to 
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FIG. 3: Radial couplings among the (4-7)"*^ Ej" states in the 
neighborhood of the n = 3 curve crossing. The notation {i,j)g 
label the radial coupling between the i-th and j-th state. 



the n — S curve crossing and thus no direct comparison 
of the radial couplings can be made. However, in the 
ungerade symmetry we can compare the energy splitting 
of the 3 and 6^S+ states at 35.9 uq. Here, we obtain 
the value 3.8 x IQ-^ Ei, which is in very good agreement 
with the value 3.9 x 10~* E\i calculated from the results 
ofRef. 36]. 



C. Adiabatic to diabatic transformation 

In order to obtain the adiabatic to diabatic transforma- 
tion (ATDT) matrix, we have numerically solved Eq. (fT5)) 
using the radial couplings described in the preceding sec- 
tion as input data. The solution has been obtained with a 
matrix version of the Runge-Kutta Fehlberg method [i^ . 
The boundary condition T = 1 has been imposed at 
i? = 50 Oq and the integration has been performed in- 
wards. As hinted above, several of the non-adiabatic 
couplings are non-vanishing as R tends to infinity. In 
practice, our choice of boundary condition corresponds 
to setting these couplings to zero beyond 50 ao. Since 
this kind of assumption has to be made in the calcula- 
tion of the scattering matrix anyway, no additional error 
is introduced due to the diabatization procedure itself. 
In the present case, the transformation to a strictly dia- 
batic basis is merely for computational reasons and is not 
discussed further. It is worth noting, however, that the 
diabatic potential energy curves we obtain here are far 
from the intuitive ones conventionally used in connection 
with the H2 system, see Figure SI 



IV. NUMERICAL PROCEDURES 

The logarithmic derivative ([55)) has been calculated by 
numerically solving the matrix Riccati equation (|36p from 
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FIG. 4: Strictly diabatic potential energy curves obtained 
from a transformation of the (1-7)^E^ adiabatic states. The 
diabatic potential energy curve associated with the H(l) + 
H(l) limit is not shown in the figure. 
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FIG. 5: Calculated cross sections for scattering into the 
11(1) + H(n), n = 1,2,3 final states. The present results 
are compared with the close-coupling results of Fussen and 
Kubach @]. 



0.54 ao out to 50 ao- For this purpose we have used an 
algorithm developed by Johnson [27|, [IMl with the grid 
size set to 5 x 10^'^ ao. Knowledge of the logarithmic 
derivative in the final integration point has allowed us 
to calculate the partial cross sections ([32]) for scattering 
within the gerade and ungerade inversion symmetries. 
The fully symmetrized integral cross sections have been 
obtained by combining the partial cross sections accord- 
ing to Eq. ((45|) . To determine at which point the series 
in Eq. (|45p could be truncated, a simple convergence cri- 
teria was introduced. This was set to terminate the sum- 
mation if the ratios of the partial cross sections and the 
accumulated integral cross sections remained less than 
10~^ for 25 terms in succession. Over the energy range 
considered here, this led to the inclusion of approximately 
250 to 3500 partial waves. 



V. RESULTS AND DISCUSSION 

A. State dependent cross sections 

The calculated cross sections for scattering into the 
H(l)-|-H(n), n — 1, 2, 3 final states are shown in Fig.O It 
is clear that the n = 3 process dominates the others over 
almost the entire energy range, which is reasonable con- 
sidering the favorable location of the n = 3 curve crossing 
at around 36 oq. In contrast, the n = 2 transitions oc- 
cur much further in and are more easily suppressed by 
the centrifugal barrier. Only for collision energies above 
5 eV, where the n — 2 cross section exhibits a mini- 
mum, is this process able to compete with neutralization 
into the n = 3 final states. This phenomena was first 
noted by Bates and Lewis [3] and has been confirmed, at 
least on a qualitative level, by almost every study since 
then 0, S @| ■ Neutralization into the n — 1 final state 



is insignificant at all energies and is included only for 
completeness. As the collision energy approaches zero, 
all cross sections gain the characteristic dependence 
that can be expected for reactions which are governed by 
the Coulomb interaction. Here, the n = 3 cross section is 
approximately a factor of 50 larger than the n = 2 cross 
section. 

Some weaker oscillations can be observed over a large 
part of the energy range. The reason for this structure 
is not completely clear, but the most likely explanation 
is in terms of Stiickelberg oscillations, i.e., in terms of 
quantum interference arising because there are several 
competing routes through the potential landscape lead- 
ing to the same asymptotic limit. 

The close-coupling results of Fussen and Kubach Q 
have also been plotted in Fig. [5] for comparison. The 
agreement between their rt = 3 cross section and the 
present one is good. The difference is about 20 percent 
at 1 eV and increases slightly with the energy. Com- 
paring the n = 2 cross sections, the two calculations are 
seen to agree well above the minimum at 5 eV, whereas 
at lower energies the differences are more pronounced. A 
rationale for this discrepancy may be found by consider- 
ing how the inversion symmetry of the electronic states 
is approached. For this purpose we show in Fig. [6] a com- 
parison of our fully symmetrized cross sections (same as 
in Fig. [5|) together with the pure gerade and ungerade 
cross sections calculated in the conventional way. We 
first note that the n = 3 cross section is fairly insensi- 
tive to the choice of inversion symmetry over the whole 
energy range. Consequently, the effect of averaging over 
these symmetries is expected to be small. A similar ob- 
servation can be made regarding the n = 2 cross section 
above 5 eV. However, at lower energies the gerade and 
ungerade results differ strongly. This indicates not only 
that our way of combining the gerade and ungerade cross 
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FIG. 6: Fully symmetrized vs. pure gerade and ungerade 
cross sections for scattering into the H(l) + H(n), n = 1, 2, 3 
final states. 

sections will affect the final result, but also that any two 
approaches that treat the inversion symmetry differently 
may lead to different outcomes. In the one electron model 
of Fussen and Kubach, this inversion symmetry is actu- 
ally broken and as a consequence their results are likely 
to differ from ours. It should be emphasized though, that 
at the energies where this happens the contribution from 
the n = 2 channels to the overall reaction is small. It 
should also be noted that the effects of treating the nu- 
clei as indistinguishable are very small, and it is equally 
good to use the simpler formula (j46|) in favor of (|45|) to 
weight the gerade and ungerade partial cross sections. 

B. Total neutralization cross section 

In Fig. [7] we show the calculated total neutralization 
cross section, i.e., the sum of the n = 1, 2, 3 cross sections. 
As could be expected from the previous subsection, this 
cross section shows a rather smooth behavior over the en- 
tire energy range. At low collision energies, the curve ap- 
proximately resembles the behavior discussed above. 
The cross section then continues to fall off until it reaches 
a broad minimum at around 20 eV. This minimum obvi- 
ously results from the nearly constant n = 3 cross section 
in combination with the rapid increase in the formation 
of n = 2 final products. 

To examine the influence of the non-vanishing asymp- 
totic couplings upon the cross section, the point at which 
the logarithmic derivative is evaluated has been varied 
over the range 50 to 80 ag. The boundary condition on 
the ATDT matrix has been varied accordingly. At the 
collision energies 0.001, 1 and 100 eV, the relative varia- 
tion in the cross section was seen to be 4 x 10^**, 7 x 10^^ 
and 1 X 10"'^, respectively. We conclude that the effects of 
the non-vanishing asymptotic couplings over the energy 
range considered here are negligible. This is consistent 
with the conclusions drawn by Borondo et al. psj . 



Additional tests have been performed in order to esti- 
mate the relevance of the radial couplings at small inter- 
nuclear distances (as discussed in Sec. IIII Ap . This has 
been done by setting all radial couplings up to a certain 
distance Rq to zero before the diabatization procedure is 
carried out. Choosing i?o = 6 ao reduces the total neu- 
tralization cross section at 0.001, 1 and 100 eV collision 
energy by 8.6, 9.9 and 5.9 percent, respectively. These 
variations in the cross section can be considered small 
(although not insignificant) and may be taken as an up- 
per limit for the influence of the higher lying ^S+^j adia- 
batic states not included in the present study, which are 
coupled to the lower lying states mainly at small inter- 
nuclear distances. A more detailed study of these effects 
is currently being pursued. 

In Fig. [7] we also show the results of earlier studies 
relevant for the mutual neutralization of and H~ at 
low collision energies. Up to the present date, no mea- 
surements of this cross section have been reported for 
collision energies below 0.15 eV. In the region 0.15 to 
3 eV, the only experimental results that are available are 
those of Moseley et al. [l3|. A comparison shows that 
these are a factor of two to three larger than the present 
data. At energies above a few eV, experiments have been 
conducted by Szucs et al. JTj and Peart and Hayton [l^] , 
among others. In the region 3 to 10 eV, their results point 
towards a cross section that is slightly lower than ours, 
which is most evident for the flrst few data points of Peart 
and Hayton, whereas for collision energies above 10 eV 
all three studies agree favorably. Comparing with earlier 
theoretical studies, our results match within 30 percent 
the close-coupling calculations of Fussen and Kubach 
and within 40 percent the Landau-Zener calculations of 
Bates and Lewis Our results also more or less over- 
lap with the Landau-Zener results obtained by Eerden et 
al. 0. 

It may be worthwhile to review the accumulated data 
on the mutual neutralization of and H" at low col- 
lision energies (below a few cV). There are now several 
calculations of this cross section that are mutually con- 
sistent within a few tens of percent. These employ such 
diverse methods as semi-classical Landau-Zener theory as 
well as model Hamiltonian and ab initio close-coupling 
schemes. In the upper part of the energy region, the 
theoretical results are in more or less agreement with 
several of the available experiments. One of these might 
even suggest a cross section that is slightly lower than 
predicted by theory. It is a fact, however, that the only 
measurement that has been published for collision ener- 
gies below 3 eV shows a cross section that is consider- 
ably higher than others reported. To clarify this situ- 
ation, new experimental efforts are obviously required. 
Recently such measurements have been undertaken by 
the group of Urbain [4§| at the Universite Catholique de 
Louvain, Belgium, using a merged beam apparatus in the 
energy range 0.001 to 0.2 eV. Although uncertainties of 
the order of 50 percent still exist, preliminary results from 
this experiment seem to support the neutralization cross 
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FIG. 7: Total cross section for the mutual neutralization of and H~. The present results are compared with the measure- 
ments by Moseley et al. Szucs et al. and Peart and Hayton as well as the calculations by Bates and Lewis [J, 
Fussen and Kubach P], and Eerden et aI7[a|. (a) Cross sections plotted over the full energy range 0.001 to 100 eV. (b) 
Enlargement of the energy range 1 to 100 eV. 



section obtained from the various theoretical treatments. 



VI. SUMMARY 



C. Rate coefficient 

At low colhsion energies, the total neutralization cross 
section may be parametrized in terms of the relative ve- 
locity V as |10| 



a{v) = Av-"^ + Bv~^ + C + Dv 



(47) 



To determine the unknown coefficients of this expression, 
we have performed a least square fit to our calculated 
cross section in the region 0.001 to 4 eV. The result is 



A 
B 
C 
D 



4.77 X 10^2 cm'^s"^ , 
-1.73 X 10"^ cm^s"^ 
1.22 X IQ-'^'^ cm^ , 
-1.57 X 10"^^ cm s , 



corresponding to a maximum error of six percent in the 
parametrized cross section. By integrating Eq. (|47p over 
a Maxwellian velocity distribution, the associated rate 
coefficient a{T) is obtained: 



a{T) 



A 



~2C 



nkT J 



B 



TT/i 



3DkT 



(48) 



where T is the ion temperature and k the Boltzmann con- 
stant. We estimate the above results to be valid over the 
temperature range 10 to 10,000 K. For comparison, the 
rate coefficient we obtain at 300 K is 1.73 x 10~^ cm'^s"^. 



In this paper we have considered the mutual neutral- 
ization of II~''-|-II~ into H(l)-|-II(n), n = 1, 2, 3 at low col- 
lision energies. To arrive at the present results, we have 
used a molecular close-coupling approach with all degrees 
of freedom treated quantum mechanically. Adiabatic po- 
tential energy curves and non-adiabatic radial couplings 
have been calculated at the FCI level of theory employing 
a large Gaussian basis set. These quantities are in good 
agreement with those obtained from more sophisticated 
electronic structure methods [M, HU, i32, 33, 34, 35, 3(|]. 
Using a strictly diabatic representation of the potential 
energy curves and coupling matrix elements, we have cal- 
culated state dependent and total neutralization cross 
sections taking into account the identity of the nuclei. 

The present results conform to the conventional view 
that for collision energies up to a few eV, almost all neu- 
tral products go into the n — S final states. Furthermore, 
a proper weighting of the gerade and ungerade cross sec- 
tions turns out to be important only for the n ~ 2 forma- 
tion at low collision energies, where these cross sections 
anyway are small. In the collision energy region below a 
few eV, our total neutralization cross section is in good 
agreement with previous theoretical studies [1, d, Q , but 
is a factor of two to three lower than that measured by 
Moseley et al 
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(A.7) 



APPENDIX 

Here, we provide analytical expressions for the solution 
matrices a;, b;, a; and /3; which are used in the defini- 
tions of the reactance and scattering matrices of Sec. lII Cl 
Open and closed channels are distinguished depending on 
whether the total energy is greater or less than the cor- 
responding threshold energy. The channels are labeled 
short-range (covalent) if 



R>p 



and Coulombic if 



" ^ ' R>p * R 



(A.l) 



(A.2) 



where Qi Q2 are the net charges situated on each 
of the nuclei. The asymptotic wave number ki for open 
channels has been defined in Eq. ([23]) . For closed chan- 
nels, 



,1/2 



(A.3) 



It is further convenient to introduce the scaled coordinate 
— kiR and the Coulomb parameter 



Vi 



MQ1Q2 



(A.4) 



We look for pairs of linearly independent solutions to 
the radial Schrodinger equation ([27|) in the region R> p. 
These are chosen such that an_i{R) and i{R) behave as 
regular and irregular solutions, and aii^i(R) and Pa^i^R) 
as incoming and outgoing wave solutions, respectively. 
For open short-range channels, two linearly independent 
solutions are the Riccati-Bessel functions of the first and 
second kind, Si{£^i) and Ci(^i), respectively We 
adopt the definitions [1^ 



aij,i 



ifcj ^ Si , 



(A.5) 



-1/2 ^ 



iSi + Cl 



(A., 



For open Coulomb channels, the radial Schrodinger equa- 
tion is solved by the regular and irregular Coulomb func- 
tions, Fi{r]i,S^i) and Gi{r]i,^i), respectively Accord- 
ingly, 



biii = Sijk^ Gi , 



(A.9) 
(A.IO) 



a,,,i=5,jk-^'^ (^^iFi + Gi^ , (A.ll) 



-1/2 ^ 



iFi + Gi 



(A.12) 



Finally, in the case of closed short-range channels, pos- 
sible solutions are £,iii{£,i) and ^iki{^i), where ii{^i) and 
ki{£,i) are the modified spherical Bessel functions of the 
first and second kind, respectively [23 | (first and third 
kind in the cited reference). Appropriate definitions 
are [H 



a,j,i - 6,j2'/\-'/''^,ii , 



= S,,2'/^7r-'/^^,ki , 



(A.13) 



(A.14) 



, = (5y(2^)i/2^,fc-i r ^ (-l)V-ifcjl , (A.15) 



(A.16) 



S,,k-'^'Ci 



Closed Coulomb channels play no role in the present 
(A. 6) study and are therefore not considered. 
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